{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 有限元算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 预备知识"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定向量函数 $\\mathbf F(x)$, 其定义域为 $\\Omega\\in\\mathbb R^n$, $\\mathbf n$ 是 $\\Omega$ 边界 $\\partial \\Omega$ 上的单位外法线向量.\n",
    "\n",
    "$$\n",
    "\\int_{\\Omega} \\nabla\\cdot\\mathbf F~ \\mathrm d x = \\int_{\\partial \\Omega}\\mathbf  F\\cdot\\mathbf n ~\\mathrm d s\n",
    "$$\n",
    "\n",
    "这就是**散度定理**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\int_{\\Omega} \\nabla\\cdot(v\\nabla u)~\\mathrm d x &= \\int_{\\partial\\Omega} v\\nabla u\\cdot\\mathbf n~\\mathrm d s\\\\\n",
    "\\int_{\\Omega} v\\Delta u~\\mathrm d x + \\int_{\\Omega}\\nabla u\\cdot\\nabla v~\\mathrm d x &= \\int_{\\partial\\Omega} v\\nabla u\\cdot\\mathbf n~\\mathrm d s\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\int_\\Omega v_x \\mathrm d x = \\int_\\Omega \\nabla\\cdot \\begin{pmatrix}\n",
    "v\\\\0\n",
    "\\end{pmatrix} \\mathrm d x = \n",
    "\\int_{\\partial \\Omega} vn_x \\mathrm d s\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用散度定理： \n",
    "* 求多边形的面积 \n",
    "* 求多面体的体积"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "考虑如下D氏边界条件的Poisson方程:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "-\\Delta u(\\mathbf x) &=& f(\\mathbf x),\\text{ on } \\Omega.\\label{eq:P}\\\\\n",
    "u|_{\\partial\\Omega} &=& 0.\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Galerkin方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有限元方法是一种基于 PDE (partial differential equations) 的变分形式 (variational formulation) 求解PDE近似解的方法."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "引入函数空间 $H^1(\\Omega)$, 对于任意 $v \\in H^1(\\Omega)$, $v$和它的一阶导数都在 $\\Omega$ 上 $L^2$ 可积. 这里的 $H^1(\\Omega)$ 是一个无限维的空间."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另外, 引入空间 $H^1_0(\\Omega) := \\{v\\in H^1(\\Omega), v|_{\\partial\\Omega} = 0\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于任意的 $v\\in H^1_0(\\Omega)$, 同乘以方程 \\eqref{eq:P} 的两端, 然后做分部积分可得: \n",
    "\n",
    "\\begin{equation}\\label{eq:wg}\n",
    "\\int_{\\Omega}\\nabla u\\cdot\\nabla v\\mathrm{d}x = \\int_{\\Omega}fv\\mathrm{d}x,\\quad\\forall v \\in H^1_0(\\Omega).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "原问题就转化为: 求解 $u\\in H_0^1(\\Omega)$, 满足\n",
    "\n",
    "\\begin{equation}\\label{eq:W}\n",
    "a(u,v) = <f,v> \\text{ for all }v\\in H_0^1(\\Omega).\n",
    "\\end{equation}\n",
    "\n",
    "其中\n",
    "\n",
    "$$\n",
    "a(u,v) = \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\mathrm{d}x,\\quad <f,v> =  \\int_{\\Omega}fv\\mathrm{d}x.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面我们考虑所谓 Galerkin 方法来求 \\eqref{eq:W} 的逼近解. 上面 $H_0^1(\\Omega)$ 是一个无限维的空间,\n",
    "为了把无限维的问题转化为有限维的问题, 引入 $H_0^1(\\Omega)$ 的一个有限维的子空间 $V$, 比如\n",
    "$V=\\mathrm{span}\\{\\phi_1,\\phi_2,\\ldots,\\phi_N\\}$. 对任何 $v \\in V$, 它都有唯一的表示\n",
    "\n",
    "$$\n",
    "v = \\sum\\limits_{i=1}^N v_i\\phi_i.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看出空间 $V$ 和 $N$ 维向量空间 $\\mathbb{R}^N$ 是同构的, 即\n",
    "\n",
    "$$\n",
    "v = \\sum\\limits_{i=1}^N v_i\\phi_i\\leftrightarrow\\mathbf{v} =\n",
    "\\begin{pmatrix}\n",
    "v_1 \\\\ v_2 \\\\ \\vdots \\\\ v_N\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "其中列向量 $\\mathbf{v}$ 是 $v$ 在基 $\\{\\phi_i\\}_{i=1}^N$ 的坐标. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面可以构造一个离散的问题: 求 $ \\tilde u = \\sum_{i=1}^{N}u_i \\phi_i \\in V$, 其对应的向量为 $\\mathbf u$, 满足\n",
    "\n",
    "\\begin{equation}\\label{eq:d}\n",
    "a(\\tilde u, v) = <f, v>,\\quad\\forall~v\\in V.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "方程 \\eqref{eq:d} 中仍然包含有无穷多个方程. 但 $V$ 是一个有限维空间, 本质上 $\\tilde u= \\sum_{i=1}^{N}u_i \\phi_i$ 只需要满足下面 $N$ 方程即可\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "a(\\tilde u, \\phi_1) = <f, \\phi_1> \\\\\n",
    "a(\\tilde u, \\phi_2) = <f, \\phi_2> \\\\\n",
    "\\vdots \\\\\n",
    "a(\\tilde u, \\phi_N) = <f, \\phi_N> \n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "即\n",
    "\\begin{cases}\n",
    "a(\\sum_{i=1}^{N}u_i \\phi_i, \\phi_1) = <f, \\phi_1> \\\\\n",
    "a(\\sum_{i=1}^{N}u_i \\phi_i, \\phi_2) = <f, \\phi_2> \\\\\n",
    "\\vdots \\\\\n",
    "a(\\sum_{i=1}^{N}u_i \\phi_i, \\phi_N) = <f, \\phi_N> \n",
    "\\end{cases}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面的 $N$ 方程可以改写为下面的形式:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "a(\\phi_1, \\phi_1) & a(\\phi_2, \\phi_1) & \\cdots & a(\\phi_N, \\phi_1) \\\\\n",
    "a(\\phi_1, \\phi_2) & a(\\phi_2, \\phi_2) & \\cdots & a(\\phi_N, \\phi_2) \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "a(\\phi_1, \\phi_N) & a(\\phi_2, \\phi_N) & \\cdots & a(\\phi_N, \\phi_N) \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_N\n",
    "\\end{pmatrix}\n",
    "= \n",
    "\\begin{pmatrix}\n",
    "<f, \\phi_1> \\\\ <f, \\phi_2> \\\\ \\vdots \\\\ <f, \\phi_N> \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "引入**刚度矩阵**(stiff matrix).\n",
    "$$\n",
    "\\mathbf{A}=(a_{ij})_{N\\times N}  \n",
    "$$\n",
    "其中 $a_{ij}=a(\\phi_i,\\phi_j)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "和**载荷矢量**(load vector) \n",
    "$$\n",
    "\\mathbf{f} = \\begin{pmatrix}\n",
    "f_1\\\\ f_2 \\\\ \\ldots \\\\f_N\n",
    "\\end{pmatrix} \n",
    "$$ \n",
    "\n",
    "其中 $ f_i=<f,\\phi_i>$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可得到如下 $N$ 阶线性方程组:\n",
    "\n",
    "$$\n",
    "\\mathbf{Au} = \\mathbf{f}.\n",
    "$$\n",
    "\n",
    "求解可得原问题的逼近解:\n",
    "\n",
    "$$\n",
    "\\tilde u = \\sum\\limits_{i=1}^N u_i\\phi_i.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 有限元方法 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 符号说明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| 符号 | 意义|\n",
    "|:------:| :----|\n",
    "| $\\Omega$ | 求解区域 |\n",
    "| $\\mathcal{T}$ | $\\Omega$ 上的三角形网格 |\n",
    "| $N$ | $\\mathcal{T}$ 的节点个数 |\n",
    "| $NC$ | $\\mathcal{T}$ 的单元个数 | \n",
    "| $x_i\\in\\mathbb{R}^2,i=1,\\ldots,N$ | 网格节点 |\n",
    "| $\\tau := (x_i,x_j,x_k)$ | $\\mathcal T$ 中由顶点 $(x_i,x_j,x_k)$ 构成三角形单元, 其中顶点按逆时针排序 |\n",
    "| $e_{ij}$ | $\\mathcal T$ 中以 $\\mathbf x_i$ 和 $\\mathbf x_j$ 为端点的一条边 |\n",
    "| $\\tau_{ij}$ | 表示边 $e_{ij}$ 从 $\\mathbf x_i$ 看向 $\\mathbf x_j$ 左手边的单元 |\n",
    "| $(l_i, l_j, l_k)$ | 顶点 $(x_i,x_j,x_k)$ 对应的三边边长 |\n",
    "| $I=\\begin{pmatrix} 1 & 0\\\\ 0& 1 \\end{pmatrix}$ | 单位矩阵 |\n",
    "| $|\\tau_m|$ |  $\\tau_m$ 的面积 |\n",
    "| $\\omega_i$ | $\\mathcal T$ 中所有以 $x_i$ 为顶点的三角形单元集合 |\n",
    "| $\\xi_{ij}$ | 边 $e_{ij}$ 相邻三角形单元的集合, 如果有两个相邻三角形单元 $\\tau_{ij}$ 和 $\\tau_{ji}$, 则为**内部边**, 只有一个相邻单元 $\\tau_{ij}$, 则为**边界边** |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 重心坐标"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定三角形单元 $\\tau$, 其三个顶点 $\\mathbf{x}_i :=(x_i,y_i)$, $\\mathbf{x}_j :=(x_j,y_j)$ 和 $\\mathbf{x}_k :=(x_k,y_k)$ 逆时针排列, 且不在同一条直线上, 那么向量 $\\vec{\\mathbf{x}_i\\mathbf{x}_j}$ 和 $\\vec{\\mathbf{x}_i\\mathbf{x}_k}$ 是线性无关的. 这等价于矩阵\n",
    "\n",
    "$$\n",
    "A = \n",
    "\\begin{pmatrix}\n",
    "x_i & x_j & x_k \\\\\n",
    "y_i & y_j & y_k \\\\\n",
    "1   & 1   & 1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "非奇异. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "任给一点 $\\mathbf{x}:=(x,y)\\in\\tau_m$, 求解下面的线性方程组\n",
    "\n",
    "$$\n",
    "A \n",
    "\\begin{pmatrix}\n",
    "\\lambda_i \\\\\n",
    "\\lambda_j\\\\\n",
    "\\lambda_k  \n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "x_i & x_j & x_k \\\\\n",
    "y_i & y_j & y_k \\\\\n",
    "1   & 1   & 1 \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\lambda_i \\\\\n",
    "\\lambda_j\\\\\n",
    "\\lambda_k  \n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "x \\\\\n",
    "y\\\\\n",
    "1  \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "可得唯一的一组解$\\lambda_i,\\lambda_j,\\lambda_k$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此对任意二维点 $\\mathbf{x}\\in\\tau$, 有\n",
    "\n",
    "$$\n",
    "\\mathbf{x}=\\lambda_i(\\mathbf{x})\\mathbf{x}_i + \\lambda_j(\\mathbf{x})\\mathbf{x}_j + \\lambda_k(\\mathbf{x})\\mathbf{x}_k \n",
    "\\text{ with } \\lambda_i(\\mathbf{x}) + \\lambda_j(\\mathbf{x}) + \\lambda_k(\\mathbf{x}) = 1. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\lambda_i,\\lambda_j,\\lambda_k$ 称为点 $\\mathbf{x}$ 关于点 $\\mathbf{x}_i,\\mathbf{x}_j$ 和$\\mathbf{x}_k$ 的**重心坐标**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "重心坐标有它相应的几何意义. 给定 $\\mathbf x\\in\\tau$, 把 $\\tau$ 的第 $i$ 个顶点 $\\mathbf{x}_i$ 换 $\\mathbf{x}$\n",
    "得到的三角形记为 $\\tau_i(\\mathbf{x})$, 则由克莱姆法则可得\n",
    "\n",
    "\\begin{equation}\\label{eq:bc}\n",
    "\\lambda_i = {|\\tau_i(\\mathbf{x})| \\over |\\tau|}.\n",
    "\\end{equation}\n",
    "\n",
    "其中 $|\\cdot|$ 表示三角形的面积."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知, $\\lambda_i, \\lambda_j, \\lambda_k$ 都是关于 $\\mathbf x$ 的线性函数, 且有\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\lambda_i(\\mathbf x_i) = 1,& \\lambda_i(\\mathbf x_j) = 0,& \\lambda_i(\\mathbf x_k) = 0\\\\\n",
    "\\lambda_j(\\mathbf x_i) = 0,& \\lambda_j(\\mathbf x_j) = 1,& \\lambda_j(\\mathbf x_k) = 0\\\\\n",
    "\\lambda_k(\\mathbf x_i) = 0,& \\lambda_k(\\mathbf x_j) = 0,& \\lambda_k(\\mathbf x_k) = 1\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\lambda_i, \\lambda_j, \\lambda_k$ 关于 $\\mathbf x$ 的梯度为:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla\\lambda_i = \\frac{1}{2|\\tau|}(\\mathbf x_k - \\mathbf x_j)W\\\\\n",
    "\\nabla\\lambda_j = \\frac{1}{2|\\tau|}(\\mathbf x_i - \\mathbf x_k)W\\\\\n",
    "\\nabla\\lambda_k = \\frac{1}{2|\\tau|}(\\mathbf x_j - \\mathbf x_i)W\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "其中 $W=\\begin{pmatrix}0& 1 \\\\-1 & 0 \\end{pmatrix}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 线性有限元基函数与空间"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定求解区域 $\\Omega$ 上的一个三角形网格 $\\mathcal T$, 它有 $N$ 个网格节点 $\\{\\mathbf x_i\\}_{i=1}^N$, $NC$ 个三角形单元 $\\{\\tau_m\\}_{m=1}^{N\n",
    "C}$.\n",
    "\n",
    "![三角形网格](figures/triangulation.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定网格节点 $\\mathbf x_i$, 记 $\\omega_i$ 为 $\\mathcal T$ 中所有以 $\\mathbf x_i$ 为顶点的三角形单元集合, 即\n",
    "\n",
    "$$\n",
    "\\omega_i = \\{\\tau,\\,\\mathbf x_i \\text{ is one vertex of }\\tau \\in \\mathcal T\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对每个网格节点 $\\mathbf x_i$, 可以定义函数 \n",
    "$$\n",
    "\\phi_i(\\mathbf x) =\n",
    "\\begin{cases}\n",
    "\\lambda_i(\\mathbf x),& \\mathbf x \\in \\tau_m \\in \\omega_i\\\\\n",
    "0, & \\mathbf x \\in \\tau_m \\notin \\omega_i\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "由 $\\phi_i(\\mathbf x)$ 的定义和重心坐标函数的性质可知:\n",
    "1. $\\phi_i(\\mathbf x_i)=1$,\n",
    "1. $\\phi_i(\\mathbf x)$ 限止在 $\\omega_i$ 中的每个单元 $\\tau$ 上, 为 $\\mathbf x_i$ 对应的重心坐标函数.\n",
    "1. $\\phi_i(\\mathbf x)$ 在 $\\omega_i$ 以外的单元上函数值为 0.\n",
    "1. $\\text{supp}(\\phi_i)=\\omega_i$.\n",
    "\n",
    "因此, 我们可以说 $\\phi_i$ 定义在 $\\mathcal T$ 上的**分片线性连续函数**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "把 $\\mathcal T$ 中的每个节点函数一起可以做为一组基, 张成一个**分片线性连续函数空间** \n",
    "$$\n",
    "V = \\text{span}\\{\\phi_1, \\phi_2, \\cdots, \\phi_N\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 刚度矩阵与右端载荷的计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "刚度矩阵 $A$ 的每一项\n",
    "$$\n",
    "a_{ij} = a(\\phi_i, \\phi_j) = \\int_{\\Omega}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x\n",
    "$$\n",
    "\n",
    "由 $\\phi_i$ 的定义可知, 我们并不需要在整个 $\\Omega$ 或者说整个 $\\mathcal T$ 上来计算上面的积分, 只需要在 $\\phi_i$ 和 $\\phi_j$ 的支集的交集上计算, 即\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\int_{\\Omega}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x\\\\\n",
    "= &\\int_{\\omega_i\\cap\\omega_j}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x\\\\\n",
    "= & \n",
    "\\begin{cases}\n",
    "\\sum_{\\tau\\in\\omega_i}\\int_{\\tau}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x, & i=j \\\\\n",
    "\\int_{\\tau_{ij}}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x + \\int_{\\tau_{ji}}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x, & e_{ij}\\text{ 为内部边} \\\\\n",
    "\\int_{\\tau_{ij}}\\nabla \\phi_i\\cdot\\nabla \\phi_j\\mathrm d\\mathbf x, & e_{ij}\\text{ 为边界边} \\\\\n",
    "0, & \\mathbf x_i \\text{ 与 } \\mathbf x_j\\text{ 不相邻}\n",
    "\\end{cases}\\\\\n",
    "= & \n",
    "\\begin{cases}\n",
    "\\sum_{\\tau\\in\\omega_i}\\int_{\\tau}\\nabla \\lambda_i\\cdot\\nabla \\lambda_j\\mathrm d\\mathbf x, & i=j \\\\\n",
    "\\int_{\\tau_{ij}}\\nabla \\lambda_i\\cdot\\nabla \\lambda_j\\mathrm d\\mathbf x + \\int_{\\tau_{ji}}\\nabla \\lambda_i\\cdot\\nabla \\lambda_j\\mathrm d\\mathbf x, & e_{ij}\\text{ 为内部边} \\\\\n",
    "\\int_{\\tau_{ij}}\\nabla \\lambda_i\\cdot\\nabla \\lambda_j\\mathrm d\\mathbf x, & e_{ij}\\text{ 为边界边} \\\\\n",
    "0, & \\mathbf x_i \\text{ 与 } \\mathbf x_j\\text{ 不相邻}\n",
    "\\end{cases}\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由以上推导可知, 我们实际上只需要在每个单元 $\\tau$ 上计算出下面六个积分, 即可组装出刚度矩阵:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\int_{\\tau}\\nabla\\lambda_i\\cdot\\nabla\\lambda_i\\,\\mathrm d \\mathbf x\\\\\n",
    "\\int_{\\tau}\\nabla\\lambda_j\\cdot\\nabla\\lambda_j\\,\\mathrm d \\mathbf x\\\\\n",
    "\\int_{\\tau}\\nabla\\lambda_k\\cdot\\nabla\\lambda_k\\,\\mathrm d \\mathbf x\\\\\n",
    "\\int_{\\tau}\\nabla\\lambda_i\\cdot\\nabla\\lambda_j\\,\\mathrm d \\mathbf x\\\\\n",
    "\\int_{\\tau}\\nabla\\lambda_i\\cdot\\nabla\\lambda_k\\,\\mathrm d \\mathbf x\\\\\n",
    "\\int_{\\tau}\\nabla\\lambda_k\\cdot\\nabla\\lambda_j\\,\\mathrm d \\mathbf x\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "上面的积分是可以直接算出来."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "类似上面分解的思想, 右端载荷向量中每一个分量, 可做如下分解:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "<f,\\phi_i> &= \\int_{\\Omega}f\\phi_i\\,\\mathrm d\\mathrm x\\\\\n",
    "& = \\int_{\\omega_i}f\\phi_i\\,\\mathrm d\\mathbf x \\\\\n",
    "& = \\sum_{\\tau\\in\\omega_i}\\int_{\\tau}f\\phi_i\\,\\mathrm d\\mathbf x\\\\\n",
    "& = \\sum_{\\tau\\in\\omega_i}\\int_{\\tau}f\\lambda_i\\,\\mathrm d\\mathbf x\n",
    "\\end{aligned}\n",
    "$$\n",
    "这意味着, 我们只要在每个三角形单元 $\\tau$ 上计算下面三个积分即可, 组装出载荷向量\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\int_{\\tau}f\\lambda_i\\,\\mathrm d\\mathbf x\\\\\n",
    "\\int_{\\tau}f\\lambda_j\\,\\mathrm d\\mathbf x\\\\\n",
    "\\int_{\\tau}f\\lambda_k\\,\\mathrm d\\mathbf x\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dirichlet 边界条件处理"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "组装出刚度矩阵 $\\mathbf A$ 和载荷向量后 $\\mathbf f$ 后, 我们还不能直接求解\n",
    "$$\n",
    "\\mathbf A\\mathbf u = \\mathbf f\n",
    "$$\n",
    "还需要进一步处理边界条件. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "模型问题中, 已经知道 $u$ 在 $\\Omega$ 边界上的值, 即边界上的网格节点处的值. 此时, 解向量 $\\mathbf u$, 可以分解为两个向量\n",
    "\n",
    "$$\n",
    "\\mathbf u = \\mathbf u_{interior} + \\mathbf u_{boundary}\n",
    "$$\n",
    "\n",
    "其中\n",
    "\n",
    "$$\n",
    "\\mathbf u_{interior}[i] = \n",
    "\\begin{cases}\n",
    "\\mathbf u_i, & \\mathbf x_i \\text{ 是一个内部点}\\\\\n",
    "0, & \\mathbf x_i \\text{ 是一个边界点}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf u_{boundary}[i] = \n",
    "\\begin{cases}\n",
    "u(x_i), & \\mathbf x_i \\text{ 是一个边界点}\\\\\n",
    "0, & \\mathbf x_i \\text{ 是一个内部点}\n",
    "\\end{cases}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "则得到的线性代数系统可以做如下的变形:\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbf A\\mathbf u &= \\mathbf f\\\\\n",
    "\\mathbf A (\\mathbf u_{interior} + \\mathbf u_{boundary}) &= \\mathbf f\\\\\n",
    "\\mathbf A \\mathbf u_{interior} &= \\mathbf f - \\mathbf A \\mathbf u_{boundary}\\\\\n",
    "\\mathbf A \\mathbf u_{interior} &= \\mathbf b\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "最后一个方程, 取 $ \\mathbf b = \\mathbf f - \\mathbf A \\mathbf u_{boundary}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于只需要求内部自由度的值, 上面最后一个方程和右端需要修改一下, 如果 $\\mathbf x_i$ 是边界点, 则\n",
    "\n",
    "1. $\\mathbf A$ 的第 $i$ 个主对角元素设为 1, 第 $i$ 行和第 $i$ 列的其它元素都设为 0, 修改后的矩阵记为 $\\bar{\\mathbf A}$.\n",
    "1. $\\mathbf b$ 的第 $i$ 个分量设为 $u(x_i)$, 修改后的右端向量记为 $\\bar{\\mathbf b}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "进而可得线性方程组\n",
    "$$\n",
    "\\bar{\\mathbf A}\\mathbf u = \\bar{\\mathbf b}\n",
    "$$\n",
    "\n",
    "注意, 如果 $x_i$ 是边界点, 上述线性方程组中的第 $i$ 个方程, 实际上就是\n",
    "\n",
    "$$\n",
    "u_i = u(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
